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Abstract An elementary introduction is presented to the study of hadron structure using 
lattice QCD. Following a brief review of relevant aspects of path integrals, the discrete lattice 
path integral is presented for gluon and quark fields and used to calculate physical observables. 
Essential aspects of instanton physics are reviewed, and it is shown how the instanton content 
is extracted from lattice gluon configurations. Finally, both comparison of results including all 
gluons with those including only instantons and the study of quark zero modes associated with 
instantons and their contributions to hadronic observables are used to show the dominant role 
of gluons in hadron structure. MIT CTP# 2701 hep-lat/9804017 

1. Introduction and Motivation 

Although a quarter of a century has passed since the experimental discovery of quarks 
and the formulation of QCD, we are only now beginning to understand the essential 
physics of the structure of light hadrons. To truly understand hadron structure, one 
must solve rather than model QCD, and the only known means to do so is the numerical 
solution of lattice field theory. But obtaining accurate numerical results for observables 
from a computer is not enough — we also need to obtain physical insight. Hence, 
our strategy is to use numerical evaluation of the QCD path integral on a lattice to 
identify the configurations that dominate the action as well as to calculate observables. 
In recent years, the algorithms and techniques of lattice QCD and the performance 
of massively parallel computers have developed to the point that we are now on the 
threshold of reliable, quantitative calculations of QCD observables. Furthermore, there 
is strong evidence from lattice calculations that the topological excitations of the gluon 
field corresponding in the semiclassical limit to instantons play a dominant role in the 
structure of light hadrons. The purpose of these lectures is to describe at an elementary 
level the basic elements of lattice QCD and how numerical solution of QCD on a lattice 
is elucidating the role of instantons in light hadrons. 

To appreciate the significance of the current lattice results, it is useful to recall the 
wide range of disparate physical pictures that have arisen from the different QCD inspired 
models introduced to model hadrons. For example, non-relativistic quark models focus 
on constituent quarks interacting via an adiabatic potential. Bag models postulate a 
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region in which relativistic current quarks are confined and interact by gluon exchange. 
Motivated by large N c arguments, Skyrme models describe the nucleon as a topological 
soliton built out of qq pairs. Finally instanton models emphasize the role of topological 
excitations of the vacuum and of the quark zero modes associated with these excitations. 

Unfortunately, phenomenology has proven inconclusive in determining which, if any, 
of these fundamentally different pictures describes the essential physics of light hadrons, 
since each of the models is rich enough that with sufficient embellishment it can be made 
to fit much of the data. Whereas perturbative QCD has proven extremely useful in ex- 
tracting quark and gluon structure functions from high energy scattering experiments, it 
is inadequate to understand their origin. Hence, it is necessary to turn to nonperturbative 
solution of QCD on the lattice. 

The physical picture that arises from this work corresponds closely to the physical 
arguments and instanton models of Shuryak and others [1, 2, 3] in which the zero modes 
associated with instantons produce localized quark states, and quark propagation pro- 
ceeds primarily by hopping between these states. Thus, QCD with light quarks is unique 
among the many-body systems with which we are familiar in the sense that the quanta 
generating the interactions cannot be subsumed into a potential but rather participate 
as essential dynamical degrees of freedom. In atoms, for example, photons play a negli- 
gible dynamical role, and to an excellent approximation may be subsumed into the static 
Coulomb potential. In nuclei, mesons play a minor dynamical role, and to a good ap- 
proximation nuclear structure maybe described in terms of two- and three-body nuclear 
forces. Indeed, experimentalists need to work very hard and pick their cases carefully to 
observe any effects of meson exchange currents. And in heavy quark systems, much of 
the physics of cc and bb bound states may be understood by subsuming the gluons into 
an adiabatic potential with Coulombic and confining behavior. The fact that nucleons 
are completely different in that gluons are crucial dynamical degrees of freedom is not 
entirely unexpected. Indeed, from perturbative QCD, we already know by the work of 
Gross and Wilczck [4] and Hoodbhoy, Ji, and Tang [5] that approximately half (16/3n/ 
to be precise, where n / is the number of active flavors and equals 5 below the top quark 
mass) of the momentum and angular momentum comes from glue in the limit of high Q 2 . 
Furthermore, experiment tells us that this behavior continues down to non-perturbative 
scales of the order of several GeV 2 . 

The discussion of lattice QCD will be based on four key, underlying ideas. The first is 
the use of path integrals. One of the great contributions of Feynman to theoretical physics 
was the formulation of quantum mechanics in terms of path integrals, which provides both 
a physical picture of quantum evolution in terms of sums of time histories and a powerful 
computational framework. For the present application, we will make use of the fact 
that the path integral eliminates the non-commuting operators of quantum mechanics or 
field theory by introducing an integral over an additional continuous variable, and thus 
effectively reduces the problem of quadrature. 

The second major idea is the introduction of Euclidean time. The basic idea is to 
write \ip) — e~P H \(j>), where e~@ H acts as a filter to project the ground state \ip) out 
of an arbitrary state \<j>) having the desired set of quantum numbers, so the continuous 
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variable in the path integral is imaginary or Euclidean time. The resulting theory has 
important connections with statistical mechanics. In the case in which one sums over a 
complete set of states and calculates the trace Tre _/3ff , one is solving field theory at finite 
temperature and /? corresponds to the physical inverse temperature. The corresponding 
path integral has the structure of classical statistical mechanics in d+1 dimensions. Many 
familiar ideas from statistical physics concerning critical behavior, order parameters, and 
Landau's theory of phase transition turn out to be useful. 

The third principal idea is lattice regularization, which replaces continuum field theory 
by a finite quantum many-body problem on a lattice. For any finite lattice spacing 
a, the maximum momentum which can arise on the lattice is p max ~ -, so that the 
lattice effectively imposes a momentum cutoff of order p max which goes to infinity as the 
lattice spacing goes to zero. One of the great accomplishments in recent years has been 
the use of renormalization group arguments and other techniques to provide convincing 
approximation to the underlying continuum theory. [7] 

The last key idea is the use of stochastic, or Monte Carlo, methods to evaluate the 
lattice path integrals. One should note at the outset that the common misnomer of Monte 
Carlo "simulations" is quite misleading. In fact, we are not simulating anything. Rather, 
we are solving an equation in the same sense as one always uses numerical analysis to 
solve equations. That is, one first selects a desired level of precision, and then using 
appropriate theorems, determines an algorithm and a number of independent samples 
which yields that precision. 

The outline of these lectures is as follows. Following this introduction, aspects of path 
integrals relevant to lattice QCD are reviewed in Section 2. The basic ideas of lattice 
QCD are presented beginning with the pure gluon sector in Section 3 and then adding 
quarks in Section 4. The role on instantons in light hadrons is discussed in the final 
section. Following an overview of instantons, it is shown how the instanton content is 
extracted from lattice gluon configurations. Results including all gluon contributions are 
compared with those including only instantons to provide one indication of the dominant 
role of instantons. Finally, direct calculation of quark zero modes, observation of quark 
localization for these modes at the locations of instantons, and demonstration that these 
modes dominate the rho and pion contributions to vector and pseudoscalar correlation 
functions provide additional indications of the role of instantons in hadron structure. 

For readers who wish to go beyond the scope of the present lecture, I recommend 
several basic references. Much of the background material is discussed in more detail 
in a text co-authored with Orland.[8] In particular, the reader is referred to Chapter 
I for treatment of coherent states and Grassmann variables, Chapter 2 for discussion 
of path integrals, and Chapter 8 for a detailed explanation of stochastic methods. A 
terse introduction to lattice gauge theory is provided by Creutz [9] and more details may 
be found in the reprint volume edited by Rebbi [10] which includes all the key articles 
through 1983 and in comprehensive texts written by Rothe [II] and by Montvay and 
Miinstcr [12]. Up-to-date reviews of recent results may be found in the proceedings of 
the yearly lattice conferences published in Nuclear Physics B Proceedings Supplements. 



4 



J.W. NEGELE 



2. Path Integrals 



2T. Feynman Path Integral. - The basic idea of the path integral is illustrated by 
considering the Feynman path integral for a single degree of freedom. The evolution 
operator e~ lHt is broken up into a large number of "time slices" separated by time 
interval e, and a complete set of states is inserted between each interval 



(2.1) 
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Then, the non-commutativity of the kinetic and potential energy operators is treated by 
the following approximation which becomes exact in the limit e — > 
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Hence the evolution operator may be expressed as the sum over all paths of the expo- 
nential of the classical action 



-%Ht\ 



i) = T>(xi, ...,x n )e 



(2.3) 



x(t)=x f 



T>(xi, ...,x n ) e* eSdassical 



,(x(t)) 



x(0) = 



The quantum mechanics of the non-commuting operators x and p has thus been rep- 
resented by an ordinary integral over an additional time variable. This result may be 
generalized to many degrees of freedom as follows 
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where a complete set of states \x\ ■ ■ ■ xn) (xi ■ ■ ■ Xn\, is inserted at each time slice. 

One important property of the path integral is that a time-ordered product is repre- 
sented as follows: 

TO(h)0(t 2 ) e- l l> dtH{t ^ = e- iH(T -^0(t 2 ) e-^-^Oih) e -«*(*i-o) 
(2.5) -> Jv(x 1 ,...,x n )e iS ^-^0(x k2 )0(x kl ) 

Hence, any path integral composed of e lS and a sequence of operators automatically 
corresponds to a time-ordered product. 
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The classical limit is obtained by including the factors of Ti which have been suppressed 
thus far and applying the stationary phase approximation 



(2.6) 
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in which case the path integral represents the sum of all quadratic fluctuations around 
the classical path. 

It is important to note that there is nothing sacred about the physical time, and 
any continuous variable may be "sliced" to treat the non-commutativity of x and p. A 
common case is the Euclidean path integral, in which real time is replaced by imaginary 
time or temperature, with the result 
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In this analytic continuation in which it — > r, the Lagrangian is effectively replaced by 
the Hamiltonian in the exponent 
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Salient properties of this Euclidean path integral are the fact that it is purely real, it has 
a well-defined measure, the Wiener measure, and it has the structure of the partition 
function of statistical mechanics with one extra dimension. 

The boundary conditions on the path integral are specified by the specific matrix 
element or elements under consideration. For example, the thermodynamic trace has the 
form 

(2.9) Tre~P H = J ' dx (a 

where 
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and thus has periodic boundary conditions. For specific matrix elements however, we 
obtain the alternative form 
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2'2. Scalar Field Theory. - Using this knowledge of the Feynman path integral, it is 
now easy to generalize to scalar field theory on a lattice. Let the continuum coordinate 
r be replaced by discrete lattice coordinates n = (ni,n 2 ,n 3 ) where the rij are integers 
and lengths will be understood to be in units of the lattice spacing a. Then one simply 
views the lattice field theory as a quantum many-body problem where the canonical 
coordinate and momentum operators x t and pi are replaced by <j>(n) and n(n) and the 
position eigenstates Xi\x) = Xi\x) are replaced by eigenstates <j>{n)\(j)) = 4>(n)\4>). On the 
spatial mesh the Hamiltonian density becomes 

J ^ r {±7r 2 (r) + ±|V#r)| 2 + V(0)} 

(2.i3) E IVw + * E w n + - ^ 2 + v W")) 1 



i=l 



where Hi denotes a displacement by one lattice site in the i th direction, J2r l 71 " 2 ^) 
corresponds to the kinetic energy . -^Pj, and the remaining terms, which we will 
denote as F[<p(n)] to avoid confusion with V[0(n)] above, correspond to a sum of one- 
and two-body potentials v(xi,Xj). Introducing time slices as before yields 



(2.i4) e -fJY.A^ 2 m+n*{n)]} = f v (0fc(n)) 



) kn \(<t> k +l(ti)-<t> k (n)f+F[<p k (n)] 



The result is a path integral defined on a four-dimensional lattice, for which we may 
introduce the obvious notation n = (no, n\, ri2, ns) where hq denotes the time label and 
rii denotes the spatial label. One observes that time slicing replaces n(n) by \<f>k+i(n) — 
(j>k(n)\ 2 = \4>(n + — 4>(n)\ 2 which has the same structure as the discrete spatial 
derivative |V</>| 2 = JZ i=1 \<p(n + /U,) — (p(n)\. Hence, a general time-ordered product 
acquires the simple form 



(2.15) TO{(j))e -pJd 3 r{^ + i(V4>f+v(4,)} ^ Jv^^Q^) 



g-SEucl. (0) 



where the Euclidean action is 



(2.16) Seuc1.(0) = E { \ E ^ n + - ^M) 2 + V } 

n I i=0 ) 

This result merits several comments. Note that Seuc1.(0) is completely symmetric in 
space and time, even though the first differences in space variables arose from a finite- 
difference approximation to the spatial derivatives whereas the time differences arose 
from the path integral time slicing. Of course, we are always free to pick different mesh 
spacings, a x and a t , in the space and time directions, respectively. Although in this 
derivation, we have gone from H to S using a discrete transfer matrix for evolution 
from one time slice to the next, it will often be useful to go backwards in the other 
direction to think of the lattice action as describing evolution of specific states under the 
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Hamiltonians H from one time slice to another in order to interpret lattice observables. 
Depending on the problem, we may be led to apply different boundary conditions in x 
and t. In the case in which all boundary conditions are periodic, the physical problem 
corresponds to finite temperature field theory in a periodic three-dimensional box and 
the shortest side of the four-dimensional box will effectively act as the temperature. 

2 '3. Coherent States. - We now need to generalize this scalar field result for general 
second quantized Fermion or Boson fields. Recall that the Feynman path integral, and 
hence the scalar field path integral, used two basic ingredients: eigenstates of x, x\x) = 
x\x), and the resolution of unity, 1 = J dx\x)(x\. The analogs of these relations for 
creation and annihilation operators are provided by Boson coherent states. 

The basic idea is seen most simply for a single creation operator a\ corresponding to 
a simple harmonic oscillator, for which 



[a.at] = i 



(2.17) | n)= 1 (at)"| ) . 



The coherent state \Z) is defined 

(2.18) l^)- e zat |o>-Efr(«Tlo> = E4il-) 

and has the following properties 

yn y(n—\) 

(2.19) a\Z) = J2 -^Mn) = n m \n 1) = Z\Z) 

(2.20) {Z \Z')=Y,{m\^r=^\n) = e z ' z ' 

mn Vm\ sfn\ 

(2.21) (Z\ : A(a\a) : \Z') = e z ' z ' A(Z* , Z') 

(2.22) j d 4£fe- z * z \Z)(Z\ = l 

The last relation is most easily demonstrated by writing the complex variable in polar 
form Z = pe 1 ^ and performing the <j) integral first. Analogous results are straightfor- 
wardly obtained [8] for a complete set of creation operators a* 

\Z) = e E a z ^\o) 
a a \Z)=Z a \Z) 
(Z\ : A{a\a) : \Z') = eS« Z « Z '°A(Z*, Z) 
(2 23) Jli^^T e~ E ° Z « Za \ Z )( Z \ = Jdn(Z)\Z)(Z\ = l 
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Proceeding as before, we obtain a path integral by time slicing 

(2.24) {Zf\e-P H \Z i ) = (Z f \e-* H J d^{Z n )\Z n ){Z n \e^ H J d^{Z n ^) ■ ■ ■ 
and the matrix element of the infinitesimal evolution operator is 

d^Z k )(Z k \e-^\Z k ^) = H dZ * k ^f ka e-^ z ^(Z k \ : e~* H ^ : +0( £ 2 )|Z fe _ 1 ) 

a 

(2.25) = n 



dZ* a dZ ka c -Y^z* ka (z ka -z lk _ 1)a )-eH(Z' ka ,z lk _ 1)a ) 
2m 



with the result 



-S(Z* ka ,Z ka ) 



(2-26) (Zf\e-^ H \Zi) = J T>(Z* ka ,Z ka ) e 

where 

(2.27) S(Z*,Z) =E4E^ - + H (Z* k , a ,Z k - l!a ) 1 

k I a J 

For Fermions with creation and annihilation operators and c Q one must take an 
additional step and introduce anticommuting Grassmann variables £, so that if c Q |£) = 
£ a |£) and C/3|£) = then we can have c a cp\£,) = £ a £[}\Q = -&£a\0 = -C0C a \Q- 

For our present purposes, one may regard this construction as a set of purely formal 
definitions. Since ^ = 0, the only allowable functions are monomials, functions are 
specified by the non- vanishing terms of their Taylor series, and the definite integral is de- 
fined by the properties / d^ a — 1 and J <i£ Q £ Q = 1. Fermion coherent states are then 
defined by 

(2.28) |0=e-5>4| ) 

and satisfy relations analogous to (2.23) and yield a path integral of the form (2.26). 
Although there are a few technical details which may be found in Ref.[8], the essential 
point is that Grassmann coherent states and path integrals have essentially the same form 
as for Bosons, except for a few crucial minus signs which do all the correct bookkeeping 
for the difference between Bosons and Fermions. 

2'4. Gaussian Integrals. - Recall the general formula for the Gaussian integral over 
complex variables 

(2.29) / TT dX * 1 dXi e-^Huxi+Jtxt+Jix; = [det^p 1 1 
J 1L 2iri 

which may be proved by changing to a basis in which H is diagonal and using J dx e~ ax2 = 
\fiTfa. An analogous result is obtained for Grassmann variables by noting that 



(2.30) 
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Hence 

(2.31) J Y[ d£d& e -«.*"««i+^+^.* = [dot H] e"'* V"' 

i 

and we see that the only difference between complex variables and Grassmann variables 
is that det.ff" appears to the power — 1 and 1, respectively. 

With this result, we are prepared to integrate out the Grassmann variables from the 
path integral. Suppose the action has the form 

(2.32) = S?M(ohij + S B (<P) 

where, for example, £*, £ might represent the Fermions tp in ^{i>— ft + m )ip + Fuv(A) 2 
and <j) represents the real Bose field A. Then 



idatM(<j>)+Sfi(<p) 



(2.33) J dCdid4>e^ M( ^ +s ^ ) = J d0e lnt 

and we are left with an integral over the real field cf) of an effective action 

(2.34) S eS (<j>) = Indet Af(^) + S B (<j>) 

In the same way, we can perform the Gaussian integrals for propagators. Consider 
first the propagator (or contraction in the language of Wick's theorem) corresponding to 
the thermodynamic average of the time-ordered product of field annihilation and creation 
operators at space-time points i = (xiU) and j — (xjtj), respectively: 

(Tipi'ipj ) = TvT^ 3 e-^ M{ ^ +SB ^ 

= j dedtd<i>te*e-t* M ^ +SB W 

(2.35) = J dcpM- 1 ((t>) ij e s <> tt ( 4 ' ) . 

The last line is obtained by differentiating Eq. (2.31) with respect to r\i and r\j which 
brings down the Grassmann variables £* and £* on the left and the inverse matrix on the 
right. The general integral with n pairs of creation and annihilation operators follows 
similarly from taking n pairs of derivatives and yields the general form of Wick's theorem: 



In 
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S Vt ■ ■ ■ SvtJVjn ■ ■ ■ S Vjl 
(2-36) = ]T(-irMr/ n . n ...Mri 
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where P denotes a permutation of the n indices. Hence, the Fermions may be integrated 
out of any physical observable when the action has the form (2.32), leaving the sum of 
all possible contractions weighted by the effective action (2.34), and we are left with an 
effective theory containing only Bosonic degrees of freedom. 

3. Lattice QCD for Gluons 

It is useful to begin the study of lattice gauge theory with the simplest possible gauge 
theory, and gradually increase the generality and complexity one step at a time. Hence, 
in this section we will completely ignore Fermions, and concentrate only on the pure 
gluon sector. This will correspond to the physical limit in which the quark mass goes 
to infinity and quarks cease to be dynamical degrees of freedom. Furthermore, we will 
begin with the simplest possible gauge group, U (1) corresponding to QED, and only after 
motivating and displaying the Wilson action for this case will we move to the non-Abelian 
SU(N) gauge theory. 

3T. U(l) Gauge Theory and the Wilson Action. - To motivate the way gauge theory 
will be formulated on a discrete space-time lattice, it is useful to recall the essential ideas 
underlying continuum gauge theory, and how the entire theory may be viewed as arising 
from the principle of gauge invariance. Therefore, let us consider a Lagrangian for a 
complex scalar field 

(3.1) L = d ll <f>*d li <l>-V(<f>*4>) 

and look for the simplest extension of the theory which is consistent with local gauge 
invariance. Note that since our final goal will be to calculate Euclidean path integrals, 
we will always write the Lagrangian and action in Euclidean form, with the result that 
9yiv = ^iiv and upper and lower Dirac indices are equivalent. Whereas L is manifestly 
invariant under the global gauge transformation <j>(x) = —e~ %a 4>{x), the derivatives in L 
yield new terms in the case of a local transformation (f>(x) = e~ ia ^(f>(x) 

(3.2) d^*8^ = - i(0„a)) 0]* (0„ - i{d„a)) <t> 

If we adopt the principle of local gauge invariance, that the theory should be independent 
of the arbitrary phase choice a(x) that various observers might choose at different points 
in space, then we may repair the theory by adding a "compensating" field A ll (x) such 
that 

(3.3) A„{x) = A^x) + -^a(x) 

If we now replace the derivative in L by the covariant derivative 

(3.4) D^(x) = (0„ + igA^x)) 4>{x) 



we observe that the transformation of A fl (x) exactly compensates for the undesircd 
derivative of a and yields an invariant Lagrangian. 
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Whereas the coupling of the new field A to <j> was determined from gauge invariancc, 
the only guiding principle for determining the action for A itself is simplicity and economy. 
Thus, we seek the simplest action for A^ involving the least number of derivatives which 
is consistent with gauge invariance and Lorentz invariance. Noting that 

(3.5) (A v + i&a) - d v (A„ + ifya) = d^A u - = 

we observe that F^ u is gauge invariant so that F^J 1 is both gauge and Lorentz invariant 
and we are thus led automatically to Maxwell's equations and the complete Lagrangian 

(3.6) L=~{d^A v -d v A^f + {D^T{D^)-V{<t>*<f>) . 

For subsequent treatment on a lattice, it is useful to note that the appropriate operator 
to compare fields at two different points x and y is the link variable 

(3.7) U(y,x)=e i f* dx " 9A » {x) 

which simply removes the arbitrary phases between the two points and yields a gauge 
invariant result. 

We now consider how to approximate this continuum theory on a space-time lattice. 
Often in numerical analysis, one may allow discrete approximations to break fundamental 
underlying symmetries. For example, when one solves the time-independent Schrodinger 
equation on a spatial mesh, one violates translational invariance. There may be small 
spurious pinning forces which reflect the fact that the energy is slightly lower when the 
solution is centered on a mesh site or centered between mesh sites, but there are no 
major qualitative errors and the quantitative errors may be strictly controlled. When 
one solves the time-dependent Schrodinger equation or time-dependent Hartree-Fock 
equation, however, one finds that it is important to enforce certain properties such as 
energy conservation and unitarity when discretizing the problem in time. In the case of 
lattice gauge theory, since by the previous argument gauge invariance plays such a crucial 
role in defining the theory, it is desirable to enforce it exactly in the lattice action. In 
contrast, as in the case of the Schrodinger equation, we will settle for an action which 
breaks Lorentz invariance, and simply insist on making the lattice spacing small enough 
that the errors are acceptably small. 

Following Wilson, we define the action in terms of directed link variables assigned to 
each of the links between sites of the space-time lattice. For U(l), we define the link 
variable from site n in the \i direction to site n + fj, as a discrete approximation to the 

integral e g X. dx A " which we denote 

(3.8) U^n) = e l8 ^ n) = e iagA ^ n +^ ^ . 

Thus 9^ (n) is a discrete approximation to g ^ x along the direction of the link, 
with A^ evaluated at the center of the link. When the direction is reversed, U n (n) — > 
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£/ n (n)t. The link variable is then a group element of U(l) and the compact variable 
9^(n) will be associated with agA fl (x) in the continuum limit. With these link variables, 
the integral over the field variables in the path integral is replaced by the invariant group 
measure for [7(f), which is ^ \\ d6. 

The fundamental building blocks of the lattice action are products of directed link 
variables taken counter-clockwise around each individual plaquette of the lattice. By 
construction, this product is gauge invariant, ensuring gauge invariance of the resulting 
action. A typical plaquette is sketched in Fig. f , fi and v are unit vectors in the horizontal 
and vertical directions, the position of the center of the plaquette is x, and we express the 
centers of each link in terms of displacements ±| in the fi or v directions. The product 
of the four group elements around the plaquette may thus be written 

□ 

_ e iag[ A n{x-% v)+A„(x+% ju)-A M (x+f v)-A„(x-% ft)] 
^ e ia 2 g{d tl A„-d,,A ll ,) 

(3.9) = e'" 2 ^" . 

Note that the discrete lattice difference operator A v [x + — A v {x — §/i) becomes the 
derivative ad^A v in the continuum limit, so that the exponent is a discrete approximation 
to the curl on the lattice and is proportional to d^A^ — d v A^ = F^ v in the continuum. 
Since each plaquette generates an approximation to , an action which corresponds to 
QED in the continuum limit may be constructed by choosing a function of U° v which 
yields F^ v plus terms that are negligible in the continuum limit. Defining the inverse 
coupling constant (3 g — l/<? 2 , where the subscript g distinguishes it from other quantities 
commonly denoted by /?, and expanding in the limit a — ► 0, the action may be written 

S(U)=(3 g J2(l-ReU D ) 
□ 

~ ^E( 1 - cos ( &2 5^)) 
~^E(^»+-) 

n{fiv} 

(3.10) -» iy d 4 xF^(x)F^(x) . 

In the second line {^iv} denotes the sum over all pairs of // and v arising from the sum 
over plaquettes and the extra factor of 1/2 in the last line accounts for the fact that each 
pair occurs twice in the double sum over repeated indices F^F^. The terms higher 
order in the lattice cutoff a vanish in the classical continuum limit and may give rise to 
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n + v 



U 3 = e i 9 aA ^( x +i' / ) 



n + p + v 



> 



JJ 4 = g^-M 1 "! M) 



U 2 = e i 9 aA "( x +^ fO 



X 



n 



U — v) 



n + fi 



Fig. 1. - An elementary plaquette of link variables. 



finite renormalization of the coupling constant in quantum field theory. The lattice gauge 
theory defined by (3.10) is in a form which may be solved directly using the Metropolis 
or heat bath methods. For readers worried that we have not systematically accounted 
for all corrections of higher order in a, a more complete derivation will be presented later 
for the general SU (N) case. 

It is useful to examine the role of Gauss' law and how the presence of external charges 
is manifested in this theory. The basic ideas are most easily sketched in the continuum 
theory. Since the Hamiltonian does not constrain the charge state of the system, we 
must project the states appearing in the path integral onto the space satisfying V • E = p 
with a specific background charge p, which may be accomplished by writing a <5-function 

in the form j Vxe $ dxdtx<y ^ '^~ p \ Remaining in temporal gauge A — and using the 
form of the path integral (2.2) in which both the coordinate x — > A and momentum 
p — > E appear, the path integral for the partition function projected onto the space with 
external source p may be written 



Equation (3.11) is an important result. Having started in temporal gauge A = 0, we 
see that enforcing Gauss' law gives rise to a projection integral over an additional field 
X which enters into the final action just like the original A field. Indeed, renaming 



(3.11) 




T>x T>A VE el dx dt t 1 - 6 ^" 3 (e 2 +b 2 )+i X (v-e- p )] 
V X VAe- 1 dxdt i^-^ x ) 2+B2 - lxp \} 



X = A so that Ai — diA$ = F oi and writing the source as a set of point charges 
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P( x ) = Hn °n${x - x n ), we obtain 

(3.12) Z= f DA tl e~I dxdtiF ^ F ^Y[e~ iqn I dtAo{Xn ' t) . 

•* n 

Thus, the Hamiltonian path integral with projection is precisely the Lagrangian path 
integral with a line of ±A fields at the positions of the fixed external ± charges. In the 
case of no external charges, we may think of the Lagrangian path integral including the 
Aq integral as the usual filter e~@ H selecting out the ground state. In the presence of 
charges, the path integral augmented by lines of A at the positions of the charges filters 
out the ground state in the presence of these sources. 

3'2. SU(N) Gauge Theory. - The generalization to non-Abelian gauge theory is 
straightforward. The link variables become group elements of SU (N) 

(3.13) U„(n) = e ia ^ XCA lW = e iagA ^ 

where the A c are Pauli matrices or Gell-Mann matrices for SU{2) or SU(3) and c is a 
color label which runs over the N 2 — 1 generators A c . The integration in the path integral 
is defined by the invariant group measure which we will denote by T>(U). 

Since derivation of the SU(N) action is one of our primary results, let us be slightly 
more careful than in the U(N) case and keep track of higher order terms in a to ascertain 
the leading error in our result. [13] As before 1 , we refer the centers of each of the links in 
Fig. 1 to the point x in the center of the plaquette, so that for SU(N) 

(3.14) U 3 = e's^f^+f") Ui = e i 9 a A,Ax-§v-) 

and the product of SU (N) group elements around an elementary plaquette is 

(3.15) U°„ = X\UiU 2 ulu\ 

□ 

With this notation, we note that U\(—a) = U\{a) and U 2 {— a) = U\(a). Hence, 
Tr U° v (-a) = Tr ulu\UiU 2 = TvU° v (a) so that TrU° v is an even function of a. Also, 
because U° v is a product of unitarity matrices, U°U°^ = 1. The continuum contri- 
bution is obtained by expanding A^(x — ^v) = A^x) — | J^^u + 9(a 2 ) as before and 

i re- 



applying the Baker-Hausdorff identity e x e g — e ( - xg+ ^^ x ^ + "') to each quantity below in 
curly brackets, with the result 



JjO^ _ | e ia fi r(A M -§a„A M + 0(a 2 )) e i9(A„ + f9 A1 A l , + C>(Q 2 ))| 

^ e iag(A tl + §d v A fl +0(a 2 )) e ig(A v -§d l ,A v +0(a 2 ))^ 



T.g{A^+A v + %{d^A v -d u A^) + i?f{A^,A,.,} + 0(a 2 )) 
x e iag(-A„-A u + ^ {d ll A u -d u A hl ) + ^[A hl .A v ]+0(a 2 )) 
(3 26) = e io-9[a{d^A„-d v A^)+iag[A^,A l ,] + 0{a 2 )} 
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Writing F^ v = d^A^—d^A^+iaglA^, A v \ and denning the next two terms in an expansion 
in powers of a as G M „ and H^, we obtain 



jjU ia 2 gF^+ia 3 G Ml ,+ia 4 H ^ 

(3.17) = 1 + ia 2 F^ + ia 3 G^ + ia 4 H^ - + 6(a 5 ) 



Unitarity of U° v implies that F, G and H are Hermitian which combined with the fact 
that Tr U° is an even function of a yields 

(3.18) Re Tr U° v = \ Tr(U° v + U°}) = Tr(l - ^g 2 F 2 u ) + 8(a 6 ) 

Defining the inverse coupling (3 g = 2N/g 2 , the SU(N) action may be written 

^) = /3 9 E( 1 -^ RcTrt/ °) 



PgC Cg 



□ 

4„2 



2N 



£ Tr(iA c F^lA"F^) + 0( fl 6 ) 



n {fj, v } 



(3.19) -^kjd^xF^ 



(x)F^(x)+0(a e ) . 



Note that summation over repeated indices is implied everywhere except where ^t^x 
denotes sums over distinct pairs, and the property TrA b A c = 2Sb c has been used in the 
second line. 

3'3. Wilson Loops and Lines. - Rephrasing the original discussion of gauge fields in 
terms of lattice variables, if there were a quark field defined on a lattice, then under 
a local gauge transformation, the field ipi at each site would be multiplied by a group 
element gi. The link variables were explicitly introduced to compensate such a gauge 
transformation, so the link variable Uij going from site i to j is multiplied by gi and gj 1 . 
Thus, the overall effect of a gauge transformation is the following: 

Uij -> giUijgj 1 
4>i -> gii>i 

(3.20) h^hg^ ■ 

In the pure gauge sector, where the only variables are link variables, it is clear from 

(3.20) that the only gauge invariant objects which can be constructed are products of 
link variables around closed paths, for which the factors of g and g^ 1 combine at each 
site. The Wilson loop is therefore defined as the trace of a closed loop of link variables 

W = Tr ft Uij 

ijec 

(3.21) =TrUijUjk---U mn U ni 
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Fig. 2. - A space-time Wilson loop defined by the chain of link variables C on a finite lattice 
(left) and a Wilson or Polyakov line defined by the chain of link variables L on a finite lattice 
(right). 



and specifies the rotation in color space that a quark would accumulate along the loop c 
from the path-ordered product P c e$ . 

To understand the physical significance of a space-time Wilson loop, it is useful to 
note that by (3.20), a quark creation operator at site j and an annihilation operator at 
site i transform under gauge transformations like a product of link variables connecting 
sites i and j 

(3.22) UikUut ■ ■ ■ U mj — > giUikUki ■ ■ ■ UmjgJ 1 . 

Thus, as far as the gluon fields are concerned, the ends of a chain of link variables are 
equivalent to an external quark-antiquark source, and the presence of such a chain of 
link variables therefore measures the response of the gluon fields to an external quark- 
antiquark source. 

Now, consider the time evolution of the system corresponding to the expectation value 
of the Wilson loop drawn in Figure 2 

JdUe-^TrUUij 
(3-23) (W) = JdUe _ 3(u ° . 

Prior to the time U, there are no color sources present, so evolution filters out the gluon 
ground state in the zero charge sector, |0) = e~ tiH \Q = 0). At time tj, the line of link 
variables between and L creates an external antiquark source at and a quark source 
at L. As discussed in connection with Eq. (3.12), the links in the time direction between 
i; and tf maintain these sources at and L. Hence, for any t between U and tf, the 
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evolution filters out the lowest gluon configuration in the presence of external quark- 
antiquark sources producing the state \ip) = e _ ( t_ *^- ff ?/;(0)^(L)|0). Finally, at time tf, 
the external sources at and L are removed by a line of links from L to 0, and the 
system is returned to the zero charge sector. Using Feynman's picturesque language of 
antiquarks corresponding to quarks propagating backwards in time, one may succinctly 
characterize the Wilson loop as measuring the response of the gluon fields to an external 
quark-like source traveling around the perimeter of the space-time loop in the direction 
of the arrows. 

Quantitatively, if tf — ti is large enough, the lowest gluon state in the presence of 
quark and antiquark sources separated by L will dominate, and (W) will be proportional 
to e~( t f~ ti ^ v ( L ^ where V(L) is the static quark-antiquark potential. Physically, this 
potential corresponds to the potential arising in heavy quark spectroscopy. Furthermore, 
at large distances, the potential in the pure gluon sector becomes linear (since the flux 
tube cannot be broken by qq pair creation) so the Wilson loop enables direct numerical 
calculation of the string tension. If the Wilson loop has I links in the time direction and 
J links in the space direction, then 

w(u) = (rr J| u} t ~ e - aIY ^ 

^ e -a 2 crIJ 
I ,J — >oo 

The exponent is thus proportional to the area for large loops, and this area behavior 
is a signature of confinement, since it arises directly from the linearly rising potential. 
Although the preceding physical argument was framed in Hamiltonian form with evo- 
lution in the time direction, it is clear that because of the symmetry of the Euclidean 
action, all space-time dimensions are equivalent and the area law reflects this symmetry. 

Since the string tension can be calculated directly from Wilson loops, it is useful to 
relate it to an experimentally measurable quantity, the slope of Regge trajectories. It 
is an empirical fact that families of meson states with a given set of internal quantum 
numbers have mass dependence on the total angular momentum J which is accurately 
described by the Regge formula 

(3.25) Mj = ^J , c/ = 0.9GeV~ 2 . 

To see how the slope a' is related to the string tension, it is useful to consider a very 
simple model in which a massless quark and antiquark are connected by a string or 
flux tube of length 2L. Since the quarks are massless, they must move at the speed of 
light, and the velocity of the segment of string a distance x from the origin is v = 
If a is defined as the energy per unit length of the flux tube in its rest frame, then the 
contributions of the element of length dx at point x to the energy and angular momentum 
are 

(3.26) dE = -yadx , dJ = jcrvxdx 



(3.24) 
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where v = (with c = 1) and 7 = (l — w 2 ) 1 ^ 2 . Hence, 

/" L adx f L a^-dx tt , 

(3.27) M=/ = 7T(tZi , J= / 9 = o^ 



and 



(3.28) M 2 = ttVL 2 = 27tctJ 
so that 

(3.29) [2to']" 1/2 = 420MeV . 

Clearly this flux tube model is a drastic oversimplification, especially for low angular 
momentum states for which the finite width of the tube, the structure of the end caps, 
and the lack of localization of the quarks could all produce large corrections. At large 
angular momentum, however, the picture is somewhat more convincing. In the context 
of this model, I will therefore take the point of view that the accurate Regge behavior 
at low angular momentum is accidental and that a is primarily determined by high 
angular momentum states. However, since the ultimate theory is undoubtedly much more 
complicated, I will not regard it as a serious problem if lattice parameters determined 
from this rough argument disagree at the 5 to 10% level relative to those determined 
from other observables such as direct calculation of hadron masses. 

Although the area law behavior of large Wilson loops is clear from (3.24), in practical 
calculations on finite lattices, there are significant corrections, including a term propor- 
tional to the perimeter arising from the self-energy of the external sources and a constant 
arising from gluon exchanges at the corners, so that 

(3.30) -InW(J, J) ~C + £>(/ + J) + a?aIJ . 



To eliminate the constant and perimeter terms, the following ratio introduced by Creutz 
[9] of Wilson loops having the same perimeter is calculated to cancel out the C and D 
terms in (3.30) 



,,„ n ,r « , ( W{I,J)W{I-\,J-\) \ 

(3 - 31) J) = ln [ w(I,J-l)W(I-l,J) ) 



2 

~ a a 



A Wilson or Polyakov line is another form of gauge invariant closed loop which can 
be placed on a periodic lattice. In this sketched in Fig. 2, the links are located 

at a fixed position in space xq and run in the time direction from the first time slice 
to the last, which by periodicity is equivalent to the first and thus renders the product 
gauge invariant. If the length of the lattice in the time direction is [3 t (were the subscript 
t distinguishes it from the inverse coupling constant /3 g = ^t), the expectation value of 
the Wilson line yields the partition function for the gluon field in the presence of a single 
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fixed quark at inverse temperature f3 t and thus specifies the free energy F quar k of a single 
quark 

/fv_ Pf^W 
V I fVUe- s ( u ) 

(3.32) = e - f3tF w' k . 

This quantity is useful as an order parameter for the deconfinement phase transition. 
Note that because the periodic lattice is a four-dimensional torus and L winds around 
the lattice once in the time direction, it is characterized by a winding number and is 
thus topologically distinct from a Wilson loop which has winding number 0. By the 
preceding argument, two lines in opposite directions, one at x — and one at x = L, 
will produce the free energy of a quark and antiquark separated by distance L, and as 
(5t — > oo this provides a means of calculating the average of the singlet and octet static 
quark-antiquark potentials. 

3'4. Strong Coupling Expansion. - One can obtain a useful physical picture of what 
happens when one evaluates the expectation value of a Wilson loop 

(3.33) (W) = Z- 1 f V(U) En^-w H u a+ U n)) trj[u 

J c 

by expanding the exponent in powers of the inverse coupling constant f3 g = ™. Since 
this expansion requires small (3 g and thus large g 2 , it is called the strong coupling expan- 
sion. Note that although it is formally analogous to the high temperature expansion in 
statistical mechanics, where f3 g would be replaced by the inverse temperature, the actual 
physical temperature of our system is specified by the length of the lattice in the time 
direction, (3 t , and is distinct from (3 g . 

The structure of the expansion is revealed by considering the integrals over group 
elements which arise in the path integral (3.33). A general discussion of integration over 
SU(N) group elements is given by Creutz [9], but for our present purposes it is sufficient 
to use the following two results, where Greek indices denote SU (N) matrix indices, not 
sites 

(3.34) J dU U af3 = 

(3.35) J dUU aSi XJ-l = IjSasSfr 

These results follow directly from the orthogonality relation for irreducible matrix repre- 
sentations of the group and are trivially verified for U(l) for which the invariant measure 
is / dU = ± \% d9 and U = e iB . 

Now consider the diagrams which result from drawing the links in the Wilson loop 
Y[ c U and some set of plaquettes U □ and Uq obtained from expanding the exponential in 
(3.33). The integral (3.34) tells us that any diagram which has a single exposed link (that 
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Fig. 3. - A 3 x 3 Wilson loop tiled with plaquettes in the strong coupling expansion. 



is, a single link between a pair of sites) anywhere on the lattice gives no contribution. 
Thus, the only non-vanishing terms in the expansion are those in which we manage to 
mate plaquettes from the exponential with the Wilson loop to eliminate all exposed links. 
The simplest way to mate two links to obtain a non-vanishing result is to place them 
between the same sites in opposite directions, which by (3.35) yields ^. Since each 
plaquette brings with it a factor of /?, the lowest order non-vanishing contribution to 
(W) is obtained by "tiling" the interior of the Wilson loop with plaquettes oriented in 
the opposite direction as sketched in Fig. 3 for a 3 x 3 loop. Note that each of the outer 
links of the original Wilson loop is protected, and the interior links protect each other 
pairwise. The leading term for an / x J Wilson loop would thus have I x J tiles, each 
contributing a factor . In addition, because of the traces in the plaquettes and Wilson 
loop in (3.33) and the 5's in Eq. (3.35) there is a factor of N for each of the (I + 1)(J+ 1) 
sites and because of the factor in (3.35) there is a factor for each of the (2IJ + I + J) 
double bonds. Hence, except for SU(2), where the counting is different because the two 
orientations of plaquettes are equivalent, the overall contribution goes as 



(3.36) 



giving the lowest order contribution to the string tension 
(3.37) a ^- a -Hn(JL) . 



Fancier tilings are also possible if one is willing to use more tiles and thus include more 
powers of (3. For example, one could place five tiles together to make a cubic box with an 
open bottom and replace one or more tiles with this box. The box could be elongated, 
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or even grown into a tube which connects back somewhere else. Alternatively, one could 
replace a plaquette oriented in one direction by (N— 1) plaquettes oriented in the opposite 
direction to obtain a non- vanishing SU(N) integral. 

The utility of this expansion is twofold. It provides a physical picture of filling in the 
Wilson loop with a gluon membrane, whose vibrations and contortions represent all the 
quantum fluctuations of the gluon field. When observed on a particular time slice, the 
cross section of this surface corresponds to a color flux tube joining the quark-antiquark 
sources. In addition, in low orders, the individual terms can be calculated explicitly and 
provide a valuable quantitative check of numerical calculations. 

3"5. Continuum Limit and Renormalization. - Pure gauge theory on a finite lattice is 
specified by two parameters: the dimensionless bare coupling constant g and the lattice 
spacing a corresponding to a momentum cutoff p max ~ f • As a is changed, the bare g 
must be changed to keep physical quantities fixed. 

In principle, the renormalization procedure on a lattice is very simple and could 
be carried out as follows. First, pick an initial value of g and calculate some set of 
dimcnsionful physical observables (Oi). These observables may be written in the form 

(3.38) (0,)=a- d *(/,( 3 )) 

where di is the dimension of the operator and /, is the dimensionless quantity calculated 
on the lattice using the Wilson action with (3 — and with all lengths expressed in 
units of the lattice spacing a. For example, we have already seen in Eq. (3.31) that the 
string tension has the form a = a~ 2 x- Then, use the physical value of one operator, 
say Oi, to determine the physical value of a corresponding to the selected g. Again, 
using the string tension example, we could define a = ^/x/420McV. With this value 
of a determined from C 1; all other observables C 2 • • • Oat are completely specified. One 
should then repeat this procedure for a sequence of successively smaller and smaller values 
of g, thereby determining the function a(g) and a sequence of values for the observables 
O2 • • • On- If the theory is correct, then each sequence of observables Oi i 7^ 1 should 
approach a limit as g — > 0, and that limit should agree with nature. 

In practice, it would be very difficult to carry out a series of calculations as described 
above to small enough g to make a convincing case. Hence, it is preferable to make 
use of our knowledge of the relation between the coupling constant and cutoff based on 
the renormalization group in the perturbative regime, and only carry out explicit lattice 
calculations down to the point at which the renormalization group behavior is clearly 
established. The foundation of the argument is the fact that the first two coefficients in 
the expansion of the renormalization group function are independent of the regular- 
ization scheme, and thus may be taken from continuum one and two loop calculations: 

(3.39) a£ = /3 g 3 + Pig 5 + ■■■ 

where (3q = y^(ll — I^f) and [}\ = j^-(102 — ^-Np). Integration of this equation 
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yields the desired relation 



(3.40) 



where is an integration constant and we have used the values of (3q and (i\ for SU (3) 
with Nf Fermions. 

The constant governing the relation between the bare coupling constant and the 
lattice cutoff can be related by one-loop continuum calculations to the constants Amom 
and A-jg^-, which govern the relation between the renormalized coupling constant and 
continuum cutoff using the momentum space subtraction procedure in Feynman gauge 
and the minimal subtraction procedure respectively, with the results [6] 



This correspondence is important for two reasons. First, the large coefficients in (3.41) 
allow us to reconcile our notion that the basic scale Aqcd is of order several hundred MeV 
with the fact that lattice measurements yield values of A^ <~ 4 — 4.6 MeV, which would 
otherwise appear astonishingly low. Second, in principle, it will provide a quantitative 
consistency test if experiments in the perturbative regime of QCD can produce sufficiently 
accurate values of Amom or A-^g. 

There is now convincing evidence that lattice calculations in the pure gauge sector 
display the correct renormalization group behavior, and thus provide accurate solutions 
of continuum QCD. Numerical evidence exists for two independent quantities, T tr , the 
temperature of the deconfinement transition, and the string tension a. The transition 
temperature on a lattice with N t time slices is given by T tI = 1/N t a(gt r ) where g tr is 
the value of the coupling for which the transition occurs. If a(g tI ) is calculated using 
the perturbative expression (3.40), then once g is small enough that the lattice theory 
coincides with the continuum theory, the quantity T tr / Kl should approach a constant. 
A number of calculations show that this is the case above (3 = 6/g 2 — 6. Similarly, 
one observes the same behavior in the string tension by plotting a/A 2 L as a function of 
(3 = 6/g 2 and observing that this ratio approaches a constant beyond (3 = 6. 

On the theoretical side, impressive progress has been made in constructing improved 
actions, which in addition to containing the leading order contribution to the action in 
Eq. 3.19, also include corrections to higher order in a (or equivalently, l/g 2 )- The most 
straightforward corrections are calculated perturbatively, but suffer from the fact that 
each additional operator in field theory has its own renormalization constant, and making 
a 10% error by calculating such a constant perturbatively is like making 10% typing 
errors in the coefficients of an allegedly high-order Runge-Kutta integration formula. 
A more effective approach is the recent use of renormalization group methods [7] to 
derive nonperturbatively corrected "perfect actions" which provide extremely impressive 
approximations to the continuum theory. 

In summary, lattice gauge theory in the pure gauge sector is quite satisfactory. There 
are no glaring conceptual or computational problems, and all the numerical evidence to 



(3.41) 



Amom = 83.5 A^ 
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date suggests that one obtains an excellent approximation to the continuum theory for 
(3 g above 6. In contrast, we will now see that full QCD including dynamical Fermions is 
more problematic at both the conceptual and computational levels. 

4. Lattice QCD with Quarks 

Significant new problems arise when one attempts to apply the ideas which are so suc- 
cessful for the stochastic evolution of path integrals for Bosons to many-Fermion systems. 
The fundamental underlying problem is the minus signs arising from antisymmetry. In 
the absence of projection onto the antisymmetric subspace, e~@ H filters out the lowest 
state of any symmetry, and favors the lowest symmetric state of energy Eq relative to the 

A f> ( fpA IT'S \ 

lowest antisymmetric state with energy Eg by the factor e~ p ^° >. If one attempts 
to project stochastically, for example by antisymmetrizing path integral Monte Carlo 
evolution at each step, the projection error is of order 1/y/N and can never overcome 
the exponential factor in cases in which Eq represents a half filled band, Fermi sea, or 
Dirac sea and Eq corresponds to a Bose condensate in the lowest state. Thus, aside from 
special cases such as one spatial dimension, the only known alternative is to write a path 
integral with Grassmann variables, introduce integrals over auxiliary fields if necessary 
to reduce it to quadratic form, and integrate out the Fermion variables as in Eq. (2.33) 
to obtain a Bosonic action containing a Fermionic determinant. Since the Grassmann 
integral has been done analytically, the projection onto the antisymmetric space is exact, 
eliminating one part of the sign problem. The resulting determinant may be positive 
or negative, so there still remains the danger of catastrophic sign cancellations in the 
stochastic evaluation of the remaining Bosonic integrals. If, however, as in our present 
case, there is an even number of Fermion species with the same action, the determinant 
appears with an even power and this final sign problem is also eliminated. 

This major detour to beat the Fermion sign problem comes at a high price. We 
started with Fermions, which in occupation number representation arc represented by a 
bunch of l's and O's, and we seek to deal with them on a digital computer which can only 
work in terms of l's and O's. Yet we must resort to calculating determinants of huge, 
non-local, nearly-singular matrices involving exceedingly large numbers of floating-point 
variables and operations. In addition, by virtue of putting the Fermions on a lattice, we 
encounter an additional unexpected difficulty associated with Fermion doubling, which 
we will discuss next. It should thus be clear at the outset, that the treatment of Fermions 
provides fertile ground for new ideas. 

4T. Naive Lattice Fermions and Doubling. - To appreciate the essential issues, con- 
sider the simplest Hermitian finite difference expression with the desired continuum limit: 



gnaivc 



a 
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J d A x 4>m^j + ^|^7p (1 + w\) (1 + ad„) V> 
- ^ I 1 + a d,j, ) 7^ (1 - igaAp) tp 



(4.1) J d 4 x^[m + ^(d^ + igA^}^ . 

Note that throughout we will use Euclidean 7-matrices satisfying 7^71, + 71, 7 M = 2<5 M „, for 
which an explicit representation is given in Ref.[9]. Although this naive action appears 
to have the desired continuum behavior and symmetries, it has an unexpected problem. 
To see this problem in its simplest form, consider the Hamiltonian corresponding to the 
action (4.1) in one space dimension for the special case of free quarks (A = 0) and zero 
mass: 



^naivc = a ^t(„) a _L {^{n + 1) _ ^{n _ 1)) 

n 

(4.2) => J dx^a-^i}) 

where a — 7071 and for convenience wc choose the representation a= [\ ^ ) . Note 



that H nmvc is Hcrmitian by virtue of the symmetric difference ip(n + 1) — ip(n — 1). We 
now transform the field operators to momentum space by writing the Fourier sum 

(4.3) V(n) = ^ E ^e ikna 

where it is understood that for a lattice with N sites and periodic boundary conditions, 
the sum over momenta in the first Brillouin zone extends over the N momenta 

The Hamiltonian is diagonal 

(4.5) H™™ = E 4a S ^^-^ k 

and thus has the eigenvalue spectrum 

with eigenfunctions 

(4.7) ^±( n )=e ikna X ± 
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^lattice 




-> x = n a 
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Fig. 4. - Fermion doubling in one dimension. The top plots compare the physical continuum 
spectrum with the spectrum of the lattice Hamiltonian. The lower plots show a half wavelength 
of the real part of the non- vanishing component of a physical wave function ^k(x) and its 
degenerate unphysical sawtooth partner ^iL-k(x). 



where \ is a two component spinor with either an upper or lower component unity and 
the other component zero. The comparison of the continuum spectrum for a massless 
Dirac particle Ek = ±fc and the lattice spectrum in the top of Fig. 4 displays the species 
doubling problem. In the region denoted by the dashed circle centered at the origin, the 
lattice spectrum (4.6) yields a good approximation to the linear physical spectrum, and 
the range of linearity increases as a — > 0. However, at the edge of the Brillouin zone, there 
is a second region in which the spectrum also goes to zero linearly, denoted by the two 
dashed semicircles. In fact, for every physical mode V&fc, there is a precisely degenerate 
unphysical mode ^n-k- Since the partition function blindly counts and weights all modes 
according to their energies, it is clear that all Fermion loops will be overcounted by a 
factor of 2 in all physical observables. Note also that since the velocity is v = the 
lattice spectrum necessarily mixes right-moving and left-moving modes. 

The origin and structure of the doubled states is simple. The degenerate partner 
to the state \&fc(n) = e lkna x is the state V^-fct™) = e m '*e~ tkna x, that is, a sawtooth 
mode in which every other lattice site has an extra factor of —1. The real part of a 
low k mode Vl/^ and its sawtooth partner ^^-k are sketched for a half wavelength in 
the lower section of Fig. 4. Note that although there arc sufficient points in the half 
wavelength of 'Ffc to yield an accurate integral with any smooth function, there is no way 
that the rapidly oscillating wave function ^n-h can represent a mode with momentum 
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near ^. Thus, we need some way to eliminate these modes so that they play no role in 
the continuum limit. The origin of the degeneracy of the physical mode with its sawtooth 
partner is the symmetric difference approximation to the derivative ip' ~ ^( n + 1 ')~^f'( n — 1) 
in the naive Hamiltonian (4.2), which clearly is impervious to the minus signs e~ m7T and 
thus yields the same magnitude for the derivatives of ipk and ip--k- The origin of this 
symmetric difference, in turn, is Hermiticity, since expressions involving only nearest 
neighbor differences like ip^n(ip(n + 1) — ip(n)) are non-Hermitian and yield complex 
eigenvalues. 

We may note in passing that a finite difference approximation to a second derivative 
strongly breaks the degeneracy between the physical and sawtooth modes. Indeed, a 
perturbation of the following form 

H' = -J2 ^ f W7o (iKn + 1) - 2%l>{n) + %l>{n - 1)) 

n 

= -«a2.fW7o( - 2 J 

n 

(4.8) — ► -a I dxi>{x)ip"{x) — ► 

a^O J a^O 

which vanishes in the continuum limit and behaves like a momentum-dependent mass 
term 

(4.9) H' = X>m(fc)*Afc 

k 

where 

2 

m(k) = - (1 — cos(ka)) 
~ ak 2 

fe^O 



(4.10) 



fc-f a 



For the physical modes in the region of the origin, the perturbation has no effect as a — > 0, 
consistent with the fact H' vanishes in the continuous limit. For the spurious sawtooth 
modes, however, the mass diverges as i. Thus, one way to remove the unphysical modes 
would be to include a term of the form (4.8) in the lattice Hamiltonian to raise their 
mass so high that they contribute negligibly to the partition function, and this is the 
idea underlying Wilson's lattice action for Fermions discussed below. 

The doubling we have discussed for simplicity in one dimension arises analogously in 
each of the four Euclidean dimensions of the naive Fermion action, Eq. (4.1), so that we 
obtain 2 4 = 16 lattice modes for each physical mode. Again, specializing to the massless 
case, the momentum space action corresponding to Eq. (4.1) is 

(4.ii) sr ve ^m^ = J2^J2 7,^^ 
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so that the inverse propagator is 

(4.12) ( T ^)- 1 =M (fc) = ^ 7 ,^^ . 

This propagator replicates the physical behavior in the region of fc M <~ fifteen times 
around points on the edge of the four-dimensional Brillouin zone at which one or more 
of the components ~ 

We are now prepared to understand both the features giving rise to the doubling 
problem and the generality of the problem. [14] Whereas the specific function sin (fc M a) in 
Eqs. (4.4a) and (4.6) is the result of using the lowest-order Hermitian difference formula 
for the derivatives in the continuum action, the most general form of the chiral symmetric, 
Hermitian action derived from discrete derivatives on a periodic lattice with the correct 
continuum limit is 

(4.13) S F = J2^kJ2^ P ^ k )^ 

k £i 

where P fJ '(k) is real for Hermiticity, P ll (k) — > for the correct continuum limit, and 

k^O 

P M (fc) is periodic under fc M — ► k^ + ^- and continuous for local discrete difference formulae 
on a lattice. Note that chiral symmetry requires the form i/^P^V, so that under a chiral 
transformation ip — > e lan " 5 ip, the two sign changes from 70 in -0 and 7^ leave the action 
invariant. Since P fJ, (k) is real, continuous and periodic in k^ with period — , it must cross 
the axis at some intermediate point, so that this general discrete action has the doubling 
observed in Eq. (4.11) in each of the four Euclidean directions yielding 15 spurious low- 
mass excitations for each physical excitation. A rigorous version of these arguments is 
known as the Nielsen-Ninomiya no-go theorem, [15] which proves using homotopy theory 
that one cannot avoid Fermion doubling in a lattice theory which is simultaneously 
Hermitian, local and chiral symmetric. 

An additional aspect of Fermion doubling is the absence of the axial anomaly. The 
axial current, which is conserved for gauge theories at the classical level but not conserved 
at the quantum level in the continuum theory, is conserved for naive lattice Fermions. 
Again, the culprits are the unphysical lattice duplicates of the physical Fermion exci- 
tations, which couple to an external axial current with the opposite chiral charge and 
effectively cancel the axial anomaly arising from the physical Fermions. 

4'2. Wilson Fermions. - One of the ways out of the no-go theorem is to give up 
chiral symmetry and, following Wilson, add a second derivative term of the form (4.8) 
to raise the mass of the unphysical sawtooth modes. In one dimension, combining the 
naive Hamiltonian (4.2) with a multiple r of the perturbation (4.8) yields the Wilson 
Hamiltonian 

~a (ip(n + 1) — ip(n — 1)) raj n il>(n + 1) — 2ip(n) + ip(n — 1)" 



n 

(4.14) = £>£[a 



2a 2i 
sin(fca) 70 (cos(fca) — l)i 

Wk 



28 



J.W. NEGELE 



Using a representation with a = 03 and —170 = 01, we obtain the energy spectrum 



(4.15) &=(*^)\\U«x{ka)-\) 
with the limits 



E -^±fc(l-ifcV + -fc 2 a 2 



(4.16) 



Thus, for fixed r, the mode for k ~ has the correct continuum limit whereas the 
sawtooth mode for fc ~ — becomes infinitely massive and decouples from the theory. 
In four Euclidean dimensions, the corresponding Wilson action is 

S w = -a 4 Y a [^(")( r - 7^(™M« + a M ) + Vi(n + a„)(r + 7^ WW] 

n H 

(4.17) + a 4 ^](m + — ^(n)ip(n) 

n 

and the propagators for the spurious modes acquire masses which diverge as - as in the 
one-dimensional case. 

The Wilson action manifestly breaks chiral symmetry for m = 0, since under the 
transformation ip — » e M75 -0, — > e l2 " 75 '0. As long as the contribution of the symme- 
try breaking term can be made arbitrarily small, its presence does not interfere with the 
physics of spontaneous symmetry breaking. In the case of a spin system, for example, 
one defines the spontaneous magnetization as the limit of the thermodynamic trace of 
the spin in the presence of an external magnetic field in the limit as the external field 
goes to zero. Heuristically, we expect the contribution of the chiral symmetry breaking 
term to vanish in the continuum limit even though the masses of the unphysical modes 
go to infinity since a J dx ^D 2 -0 — ► 0, and explicit calculation [14] shows that this is the 

a — >0 

case. Furthermore, including 15 massive Wilson Fermion partners as well as the physical 
mode yields the correct axial anomaly in the continuum limit. [14] 

The combination M = (m + — ) in Eq. (4.17) enters the action like a mass term. 
Since this mass term is not protected from renormalization by any symmetry, it must 
be renormalized by "fine tuning." Thus, by a search involving a series of calculations of 
the pion mass at fixed coupling constant g and Wilson r, we may determine a critical 
mass M CI (g,r) such that = for a chiral symmetric theory and M(g,r) such that 
TOtt = 140 MeV for the physical theory. To the extent to which it is meaningful to define 
a quark mass, one may define M quar k = M — M cr . A more conventional notation is in 
terms of the hopping parameter 



(4.18) 



1 



K = 



2Ma 2ma + 8r 
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and rescaled Fermion fields 

(4.19) * ee (Ma 4 ) 1/2 V 
for which the action has the form 

n ^ fJ. 

(4.20) +*(n + a, t )(r + 7 M )C/,I(n)*(n)]J 
and 

The fields have been scaled such that the diagonal term is now unity, and the hopping 
parameter k specifies the strength of the nearest-neighbor coupling via link variables. 
Although the Wilson parameter r is often chosen to be 1, in principle it may be optimized 
to render the errors from aripD rip in the physical modes and the contribution of the 
sawtooth modes comparable. 

Just as one can derive the conserved vector current of the continuum theory as the 
Noether current of the continuum action, one can derive a discrete version of Noether's 
theorem from the Wilson action (4.17) of the form 

(4.22) A^(n) = V^(n)-V fl (n-a fl ) = 
where the conserved vector current on the lattice is 

Vf,(n) = -\ip{n)(r - 7 /J )f7 jU (n)V>(n + a M ) 

(4.23) +^(3 + ^(r + 7^W(3) • 

Note that neither the local current ip(n)j lx ip(n) nor the point-split current defined by 
(4.23) with r = is conserved on the lattice. In contrast, since there is no chiral 
symmetry, there is no unique definition of the lattice axial current, so it is necessary to 
use perturbation theory to explicitly calculate the difference between any choice for the 
lattice axial current and the continuum axial current. 

In summary, Wilson Fcrmions provide one possible framework for solving the dou- 
bling problem on the lattice and calculations discussed subsequently will be based on this 
formulation. Other alternatives include staggered fcrmions [16], which while not avoiding 
the no-go theorem, maintain a remnant of chiral symmetry while thinning the spurious de- 
grees of freedom from 16 down to 4 in four dimensions, and Kaplan-Shamir fermions [17], 
which arc formulated on the four-dimensional boundary of a five-dimensional sphere. 



4'3. Hopping parameter expansion. - Just as the strong-coupling expansion provided 
insight into the solutions of lattice QCD in the pure gauge sector, expansion in powers 
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of the hopping parameter, K, provides analogous insights into solutions in the presence 
of Fermions. The basic idea is to expand the integral 



in powers of the hopping term xj)KUip n ±i- Integration over all ^ip then yields the sum 
of all possible closed chains of KU in which the f/'s are oriented head to tail, and these 
closed quark paths represent all the quark time histories. As before, integrating over all 
U then tiles the closed quark loops with gluons, giving rise to all the color singlet flux 
tubes connecting the quarks and antiquarks on any time slice. 

Now, consider a path integral corresponding to the propagation of a meson from 
space-time point x to y, which we write schematically, ignoring 7 matrices, as 



where M = (1 + KU). Expanding M~J-(U) = (1 + KU)~^ in powers of K generates all 
valence quark trajectories from x to y and similarly M~y(U) generates trajectories from 
y to x. 

Expansion of lndetM([7) generates all disconnected quark loops corresponding to 
excitation of quark- antiquark pairs from the Dirac sea, and expansion of S(U) tiles sur- 
faces between these valence and sea-quark trajectories. A typical configuration showing 
a valence quark-antiquark pair propagating from x to y, a quark-antiquark loop excited 
out of the fermi sea, and a minimal tiling by plaquettes is shown in Fig. 5. Cutting 
this tiling at various time slices corresponds to a quark-antiquark pair connected by a 
flux tube or two quark-antiquark pairs connected by flux tubes. From this argument, 
we see that omission of the determinant, which is very expensive computationally, yields 
the quenched or valence approximation in which quark-antiquark pairs excited from the 
sea are neglected. Typical lattices used in numerical calculation vary from 16 2 x 32 to 
32 3 x 64 sites and thus involve integration over 10 7 to 10 8 real variables. 

4'4. Correlation functions. - As in the case of other strongly interacting many-body 
systems, to understand the structure of the vacuum and light hadrons in nonperturbative 
QCD, it is instructive to study appropriately selected ground state correlation functions, 
to calculate their properties quantitatively, and to understand their behavior physically. 

Because of our subsequent interest in instantons, we will focus our attention on vac- 
uum point-to-point equal time correlation functions of hadronic currents 



discussed in detail by Shuryak [1] and recently calculated in quenched lattice QCD. [18] 
The motivation for supplementing knowledge of hadron bound state properties by these 



(4.24) 




(4.25) 




(4.26) 



R(x) = (n\TJ(x)J(0)\ST) 
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Fig. 5. - A typical time history arising in the hopping parameter expansion. The heavy lines 
connecting x and y denote valence quarks, the interior heavy line corresponds to a quark- 
antiquark pair, and the light lines denote gluons. 

correlation functions is clear if one considers the deuteron. Simply knowing the binding 
energy, rms radius, quadruple moment and other ground state properties yields very lit- 
tle information about the nucleon-nucleon interaction in each spin, isospin and angular 
momentum channel as a function of spatial separation. To understand the nuclear inter- 
action in detail, one inevitably would be led to study nucleon-nucleon scattering phase 
shifts. Although, regrettably, our experimental colleagues have been most inept in pro- 
viding us with quark- antiquark phase shifts, the same physical information is contained 
in the vacuum hadron current correlation functions R(x). As shown by Shuryak [1], 
in many channels these correlators may be determined or significantly constrained from 
experimental data using dispersion relations. Since numerical calculations on the lattice 
agree with empirical results where available, we regard the lattice results as valid solu- 
tions of QCD in all channels and thus use them to obtain information comparable to 
scattering phase shifts. 

The correlation functions we calculate in the pseudoscalar, vector, nucleon and Delta 
channels are 

R{x) = (n\TJ p (x)J p (0)\rt) 
R(x) = (filTJ^aOJ^OJIfi) 
R(x) = \Tr((n\TJ N (x)J N (0)\Q)x^ v ) 
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and 



R(x) = ±Tr 



((n\Tj£(x)j*(o)\n)x^) 



where 



.P = 



uj 5 d 




and 




e abc [u a C^U b }u C . 



As in Refs. [1] and [18], we consider the ratio of the correlation function in QCD to the 
correlation function for non-interacting massless quarks, R(x)/ Rq(x), which approaches 
one as x — > and displays a broad range of non-perturbative effects for x of the order of 
1 fm. Typical results of lattice calculations of ratios of vacuum correlation functions are 
shown in Fig. 6. 

Note that the lattice results (solid line) agree well with phcnomenological results from 
dispersion analysis of data (long dashed curves). Also, observe that the vector and pseu- 
doscalar correlation functions are strongly dominated by the rho and pion contributions 
(dotted lines) in the region of 0.5 to 1.5 fm. We will subsequently show that these rho 
and pion contributions in turn arise from the zero mode contributions associated with 
instantons. 

As discussed in Refs. [1, 18], these vacuum correlators show strong indications of in- 
stanton dominated physics. As shown by 't Hooft [19], the instanton induced interaction 
couples quarks and antiquarks of opposite chirality leading to strong attractive and re- 
pulsive forces in the pseudoscalar and scalar channels respectively and no interaction to 
leading order in the vector channel. Just this qualitative behavior is observed at short 
distance in all the channels we computed. Furthermore, as shown by the open circles 
with error bars in Fig. 6, the random instanton model of Shuryak et al. [2] reproduces 
the main features of the correlation functions at large distance as well. 

5. The Role of Instantons in Light Hadrons 

Having established the framework of lattice QCD, I will now use it as a tool to elu- 
cidate the role of instantons in light hadrons. The QCD vacuum is understood as a 
superposition of an infinite number of states of different winding number, where the 
winding number characterizes the number of times the group manifold is covered when 
one covers the physical space. Just as there is a stationary point in the action of the 
Euclidean Feynman path integral for a double well potential corresponding to the tun- 
neling between the two degenerate minima, so also there is a classical solution to the 
QCD equations in Euclidean time, known as an instanton [20], which describes tunneling 
between two vacuum states of differing winding number. The action associated with an 
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Fig. 6. - Vector (V) and Pseudoscalar (P) correlation functions are shown in the upper and 
lower panels respectively. Lattice results [18] are denoted by the solid points with error bars and 
fit by the solid curves, which may be decomposed into continuum and resonance components 
denoted by short dashed and dotted curves respectively. Phenomenological results determined 
by dispersion analysis of experimental data in Ref. [1] are shown by long dashed curves, and the 
open circles denote the results of the random instanton model of Ref. [2]. 



instanton is 

(5.1) So - 1/ A W - £ J Ar (j^)* - f . 

Note that the action density has a universal shape characterized by a size p, and that 
the action is independent of p. Furthermore, the instanton field strength is self-dual, i.e. 
F* = efivapF^p — ±F" , so that the topological change of an instanton is 

(5.2) Q = £ti J ^xF^F% - ±1 . 
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Two features of instantons are particularly relevant to light hadron physics. The first is 
the fact that although the fcrmion spectrum is identical at each minimum of the vacuum, 
quarks of opposite chirality are raised or lowered one level between adjacent minima. 
Thus, an instanton absorbs a left-handed quark of each flavor and emits a right-handed 
quark of each flavor, and an anti-instanton absorbs right-handed quarks and emits left- 
handed quarks. Omitting heavier quarks for simplicity, the resulting 't Hooft interaction 
involving the operator URU^dndLSns^ is the natural mechanism to describe otherwise 
puzzling aspects of light hadrons. It is the natural mechanism to flip the helicity of a 
valence quark and transmit this helicity to the glue and quark-antiquark pairs, thereby 
explaining the so-called "spin crisis." It also explains why the two valence u quarks in the 
proton would induce twice as many dd pairs as the uu pairs induced by the single valence 
d quark. The second feature is that each instanton gives rise to a localized zero mode of 
the Dirac operator D^ fl (p (x) = 0. Hence, considering a spectral representation of the 
quark propagator, it is natural that the propagator for the light quarks is dominated by 
these zero modes at low energy. This gives rise to a physical picture in which qq pairs 
propagate by "hopping" between localized modes associated with instantons. 

5T. Identifying instantons by cooling. - The Feynman path integral for a quantum 
mechanical problem with degenerate minima is dominated by paths that fluctuate around 
stationary solutions to the classical Euclidean action connecting these minima. [8] In the 
case of the double well potential, a typical Feynman path is composed of segments fluc- 
tuating around the left and right minima joined by segments crossing the barrier. If one 
had such a trajectory as an initial condition, one could find the nearest stationary solu- 
tion to the classical action numerically by using an iterative local relaxation algorithm. 
In this method, which has come to be known as cooling, one sequentially minimizes the 
action locally as a function of the coordinate on each time slice and iteratively approaches 
a stationary solution. In the case of the double well, the trajectory approaches straight 
lines in the two minima joined by kinks and anti-kinks crossing the barrier and the struc- 
ture of the trajectory can be characterized by the number and positions of the kinks and 
anti-kinks. 

In QCD, the corresponding classical stationary solutions to the Euclidean action for 
the gauge field connecting degenerate minima of the vacuum are instantons, and we 
apply the analogous cooling technique [21] to identify the instantons corresponding to 
each gauge field configuration. 

The results of using 25 cooling steps as a filter to extract the instanton content of a 
typical gluon configuration are shown in Fig. 7, taken from Ref. [22] using the Wilson 
action on a 16 3 x 24 lattice at 6/g 2 = 5.7. As one can see, there is no recognizable 
structure before cooling. Large, short wavelength fluctuations of the order of the lattice 
spacing dominate both the action and topological charge density. After 25 cooling steps, 
three instantons and two anti-instantons can be identified clearly. The action density 
peaks are completely correlated in position and shape with the topological charge density 
peaks for instantons and with the topological charge density valleys for anti-instantons. 
Note that both the action and topological charge densities are reduced by more than 
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two orders of magnitude so that the fluctuations removed by cooling are several orders 
of magnitude larger than the topological excitations that are retained. 




Fig. 7. - Instanton content of a typical slice of a gluon configuration at fixed x and y as a 
function of z and t. The left column shows the action density S(l, 1, z, i) before cooling (a) and 
after cooling for 25 steps (c). The right column shows the topological charge density Q(l, 1, z, t) 
before cooling (b) and after cooling for 25 steps. 



Setting the coupling constant, or equivalcntly, the lattice spacing, and quark mass by 
the nucleon and pion masses in the usual way, it turns out that the characteristic size of 
the instantons identified by cooling is 0.36 fm and the density is 1.6 fm -4 , in reasonable 
agreement with the value of 0.33 fm and 1.0 fm~ 4 in the liquid instanton model. [2] 

5'2. Comparison of results with all gluons and with only instantons. - One dramatic 
indication of the role of instantons in light hadrons is to compare observables calculated 
using all gluon contributions with those obtained using only the instantons remaining 
after cooling. Note that there are truly dramatic differences in the gluon content before 
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and after cooling. Not only has the action density decreased by two orders of magnitude, 
but also the string tension has decreased to 27% of its original value and the Coulombic 
and magnetic hypcrfinc components of the quark-quark potential are essentially zero. 
Hence, for example, the energies and wave functions of charmed and B mesons would be 
drastically changed. 
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Fig. 8. - Comparison of rho observables calculated with all gluon configurations and only 
instantons. The upper left-hand plot shows the vacuum correlator in the rho channel calculated 
with all gluons as in Fig. 6 and the upper right-hand plot shows the analogous result with only 
instantons. The lower plot shows the ground state density-density correlation function for the 
rho with all gluons (solid circles) and with only instantons (open circles). Error bars for the 
solid circles are comparable to the open circles and have been suppressed for clarity. 



As shown in Fig. 8, however, the properties of the rho meson are virtually unchanged. 
The vacuum correlation function in the rho (vector) channel and the spatial distribution 
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of the quarks in the rho ground state, given by the ground state density-density correla- 
tion function [23] (p\qjoq(x)qjoq(0)\p) , are statistically indistinguishable before and after 
cooling. Also, as shown in Ref. [22], the rho mass is unchanged within its 10% statistical 
error. In addition, the pseudoscalar, nucleon, and delta vacuum correlation functions and 
nucleon and pion density-density correlation functions are also qualitatively unchanged 
after cooling, except for the removal of the small Coulomb induced cusp at the origin of 
the pion. 

Although these cooling studies strongly indicate that instantons play an essential role 
in light quark physics, cooling has the disadvantage of modifying the instanton content of 
the original gluon configuration. It is possible to avoid the gradual shrinkage of a single 
instanton until it eventually falls through the lattice by using an improved action that is 
sufficiently scale independent. [24] However, pairs of instantons and anti-instantons will 
eventually attract each other and annihilate, thereby continually eroding the original 
distribution. Hence, it is valuable to complement these cooling calculations by studies 
of the zero modes associated with instantons, which, as we show in the next section, can 
be carried out successfully on the original uncooled gluon configurations. 

5'3. Eigenmodes of the Dirac operator. - In the continuum limit, the Dirac operator 
for Wilson fermions approaches the familiar continuum result 

Dip x = ip x - K ^2 [i r - l^Ux^ipx+t, + {r + 7(t)4- t i,)i^-(< 

(5.3) -> l[ m + i0 + 5 /L)]^ • 

In the free case, the continuum spectrum is ^ [m+i|p|] and the Wilson lattice operator 
approximates this spectrum in the physical regime and pushes the unphysical fermion 
modes to very large (real) masses. In the presence of an instanton of size p at x = 0, it 
is shown in Ref. [25] that the lattice operator produces a mode with zero imaginary part 
that approaches the continuum result 

(5.4) 1po(x)s,a=U a>a ^~ 



7T (x 2 +p 2 ) 3 / 2 

and whose mixing with other modes goes to zero as the lattice volume goes to infinity. In 
addition, instanton-anti-instanton pairs that interact sufficiently form complex conjugate 
pairs of eigenvalues that move slightly off the real axis. Thus, by observing the Dirac 
spectrum for a lattice gluon configuration containing a collection of instantons and anti- 
instantons, it is possible to identify zero modes directly in the spectrum. 

Fig. 9 shows the lowest 64 complex eigenvalues of the Dirac operator on a 16 4 un- 
quenched gluon configuration for 6/g 2 = 5.5 and k = 0.16, both before and after cooling 
(where 100 relaxation steps with a parallel algorithm are comparable to 25 cooling steps). 
The lower, cooled, plot has just the structure we expect with a number of isolated in- 
stantons with modes on the real axis and pairs of interacting instantons slightly off the 
real axis. However, even though the uncooled case shown in the upper plot also contains 
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Fig. 9. - Lowest 64 complex eigenvalues of the Wilson-Dirac operator for an unquenched gluon 
configuration both before (upper plot) and after cooling (lower plot). The scale is such that 
0.06 on the imaginary axis roughly corresponds to the lowest Matsubara frequency, 380 MeV. 



fluctuations several orders of magnitude larger than the instantons (as seen in Fig. 7) , it 
shows the same structure of isolated instantons and interacting pairs. To set the scale, 
note that if we had antiperiodic boundary conditions in time, the lowest Matsubara mode 
(ip = ij) would occur at 0.06 on the imaginary axis, so all the modes below this value 
are presumably the results of zero modes. 
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5'4. Zero mode expansion. - The Wilson-Dirac operator has the property that D = 



where A, = A* . A clear indication of the role of zero modes in light hadron observables 
is the degree to which truncation of the expansion to the zero mode zone reproduces the 
result with the complete propagator. 

Fig. 10 shows the result of truncating the vacuum correlation functions for the vector 
and pseudoscalar channels to include only low eigenmodcs.[25] On a 16 4 lattice, the full 
propagator contains 786,432 modes. The top plot of Fig. 10 shows the result of including 
the lowest 16, 32, 64, 96, and finally 128 modes. Note that the first 64 modes reproduce 
most of the strength in the rho resonance peak pointed out in Fig. 6, and by the time 
we include the first 128 modes, all the strength is accounted for. Similarly, the lower 
plot in Fig. 10 shows that the lowest 128 modes also account for the analogous pion 
contribution to the pseudoscalar vacuum correlation function. Thus, without having to 
resort to cooling, by looking directly at the contribution of the lowest eigenf unctions, we 
have shown that the zero modes associated with instantons dominate the propagation of 
rho and pi mesons in the QCD vacuum. 

5 "5. Localization. - Finally, it is interesting to ask whether the lattice zero mode 
eigenfunctions are localized on instantons. This was studied by plotting the quark den- 
sity distribution for individual eigenmodes in the x-z plane for all values of y and t, and 
comparing with analogous plots of the action density. As expected, for a cooled con- 
figuration the eigenmodes correspond to linear combinations of localized zero modes at 
each of the instantons. (Because there are no symmetries, the coefficients are much more 
complicated than the even and odd combinations in a double well or the Bloch waves in 
a periodic potential.) What is truly remarkable, however, is that the eigenfunctions of 
the uncooled configurations also exhibit localized peaks at locations at which instantons 
are identified by cooling. Thus, in spite of the fluctuations several orders of magnitude 
larger than the instanton fields themselves, the light quarks essentially average out these 
fluctuations and produce localized peaks at the topological excitations. 

6. Conclusion 

Taken as a whole, lattice calculations now provide strong evidence that instantons 
play a dominant role in quark propagation in the vacuum and in light hadron struc- 
ture. The instanton content of gluon configurations has been extracted by cooling, and 
the instanton size and density are consistent with the instanton liquid model. Vacuum 
correlation functions, ground state density-density correlation functions, and masses cal- 
culated with only instantons show striking agreement with the results obtained with all 
gluons. Zero modes associated with instantons are clearly evident in the Dirac spectrum, 



75-DT7 5 , which implies that (V'jlTsl^/'i) 
representation of the propagator 



unless Ai = A* and we may write the spectral 
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Fig. 10. - Contributions of low Dirac eigenmodes to the vector (upper graph) and pseudoscalar 
(lower graph) vacuum correlation functions. The upper graph shows the contributions of 16, 
32, 64, 96, and 128 eigenmodes compared with the full correlation function for an unquenched 
configuration with a 63 MeV valence quark mass. The lower graph compares 128 eigenmodes 
with the full correlation function for a quenched configuration with a 23 MeV quark mass. 
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and account for the rho and pi contributions to vector and pseudoscalar vacuum corre- 
lation functions. Finally, quark localization at instantons has been observed directly in 
uncooled configurations. Hence, in addition to providing a powerful tool for calculating 
a variety of important physical observables, ranging from spectroscopy and weak ma- 
trix elements to QCD thermodynamics, Lattice QCD is also beginning to teach us the 
underlying physics of hadron structure, bit by bit. 
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